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ABSTRACT 

In this paper we extend a previous work (j Vielva fc Sanz|[2009 ) where we presented a 
method based on the N-point probabihty distribution (pdf) to study the Gaussianity 
of the cosmic microwave background (CMB). We explore a local non- linear perturba- 
tive model up to third order as a general characterization of the CMB anisotropies. 
We focus our analysis in large scale anisotropies (6* > 1°). At these angular scales (the 
Sachs- Wolfe regime), the non-Gaussian description proposed in this work defaults (un- 
der certain conditions) to an approximated local form of the weak non-linear coupling 
inflationary model. In particular, the quadratic and cubic terms are governed by the 
non-linear coupling parameters fNL and gNL, respectively. The extension proposed in 
this paper allows us to directly constrain these non-linear parameters. Applying the 
proposed methodology to WMAP 5-yr data, we obtain —5.6 x 10^ < gNL < 6.4 x 10^, 
at 95% CL. This result is in agreement with previous findings obtained for equivalent 
non-Gaussian models and with different non-Gaussian estimators, although this is the 
first direct constrain on gNL from CMB data. A model selection test is performed, 
indicating that a Gaussian model (i.e. fNL = and gNL = ) is preferred to the non- 
Gaussian scenario. When comparing different non-Gaussian models, we observe that 
a pure fNL model (i.e. gNL = 0) is the most favoured case, and that a pure gNL model 
(i.e. fNL = 0) is more likely than a general non-Gaussian scenario (i.e. fNL 7^ and 
gNL 7^ 0). Finally, we have analyzed the WMAP data in two ind ependent hemispheres , 
in particular the ones defined by the dipolar pattern found by Hoftuft et al. (2009). 
We show that, whereas the gNL value is compatible between both hemisp heres, it is 
not the c ase for fNL (with a p- value « 0.04). However, if, as suggested by Hoftuft et 
ar](2009), anisotropy of the data is assumed, the distance between the likelihood dis- 
triDutions for each hemisphere is larger than expected from Gaussian and anisotropic 
simulations, not only for fNL, but also for g„Z (with a p-value of ~ 0.001 in the case 
of this latter parameter) . This result is an extra evidence for the CMB asymmetries 
previously reported in WMAP data. 

Key words: cosmology: observations - cosmology: cosmic microwave background - 
methods: data analysis - methods: statistical 



1 INTRODUCTION 

Current observations of the Cosmic Microwave Background 
(CMB) temperature and polarization fluctuations, in addi- 
tion to other astronomical data sets (e.g. [Komatsu et al. 
[2009i [Gupta et aLlpOOQ] see |Barreiro||2009f for a recent re- 
view), provide an overall picture for the origin, evolution, 
and matter and energy content of the universe, which is usu- 
ally referred to as the standard cosmological model. In this 
context, we believe the universe to be highly homogeneous 
and isotropic, in expansion, well described by a Friedmman- 



Robertson- Walker metric and with a trivial topology. The 
space geometry is very close to flat, and it is flUed with cold 
dark matter (CDM) and dark energy (in the form of a cos- 
mological constant. A), in addition to baryonic matter and 
electromagnetic radiation. Large scale structure (LSS) is as- 
sumed to be formed by the gravitational collapse of an ini- 
tially smooth distribution of adiabatic matter fluctuations, 
which were seeded by initial Gaussian quantum fluctuations 
generated in a very early inflationary stage of the universe 
evolution. 

It is interesting to mention that, besides the success of 
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current high quahty CMB data (in particular the data pro- 
vided by the WMAP satellite |Hinshaw et al.||2009[) in con- 



straining the cosmological parameters with good accuracy 
and in showing the high degree of homogeneity and isotropy 
of the Universe (as predicted by the standard inflation sce- 
nario, see for instance Liddle & Lyth 20001, it has been. 



precisely, through the analysis of this very same data, that 
the CMB community has been allowed to probe fundamen- 
tal principles and assumptions of the standard cosmological 
model. Most notably, the application of sophisticated sta- 
tistical analysis to CMB data might help us to understand 
whether the temperature fluctuations of the primordial ra- 
diation are compatible with the fundamental isotropic and 
Gaussian standard predictions from the inflationary phase. 

Indeed, the interest of the cosmology community in this 
fleld has experienced a signiflcant growth, since several anal- 
yses of the WMAP data have reported some hints for depar- 
ture from isotropy and Gaussianity of the CMB temperature 
distribution. The literature on the subject is very large, and 
is still growing, which makes really difficult to provide a 
complete and updated list of publications. We refer to our 



previous work (Vielva & Sanz 20091 for an (almost) com- 
plete review. 

Among the previously mentioned analyses, those re- 
lated to the study of non-standard inflationary models have 
attracted a larger attention. For instance, the non-linear 
coupling parameter fNL that describes the non-linear evo- 



lution of the inflationary potential (see e.g. Bartolo et al. 
2004 and references therein) has been constrained by sev- 
eral groups, from the analysis of the WMAP data: using the 
angular bispectrum (Komatsu et al.||2003 



2007; Spergel etal 



2007 Komatsu et a' 



Creminelli et al. 



.| |2009| lYadaTfc 



Wandelt,|2008 Smith et al. 


[2009 


1; applying the Minkowski 


functionals (Komatsu et al. 


2003 


Spergel et al.||2007 Gott 


et al.||2007 Hikage et al.||2008 Komatsu et al.||2009 1; using 


different statistics based on wavelets (|Mukherjee & Wang 


20041 


Cabella et al.||2005| 


:]urto et al.||2009a|b| |Pietrobon 


et al. 


2009 Rudjord et al. 


20091, and by exploring the N- 



fNL > (with a probability of around 95%, Yadav & Wan- 
delt|2008 Rudjord et al.|2009 1 , there is a general consensus 
on the WMAP compatibility with the predictions made by 
the standard inflationary scenario at least at 95% confldence 
level. The current best limits obtained from the CMB data 
are: -4 < fNL < 80 at 95% CL by 'S mith et"ar] (|2009). In ad- 
dition, very recently, promising constraints coming out from 
the analysis of LSS have been reported: —29 < fNL < 70 at 
95% CL ( [Slozar et al.pOOSt . 

The aim of this paper is to extend our previous 
work (Vielva & Sanz ] [2009| , where the full N-pdf of a non- 
Gaussian model that describes the CMB anisotropics as a 
local (pixel-by-pixel) non-linear expansion of the temper- 
ature fluctuations (up to second order) was derived. For 
this model — that, at large scales, can be considered as an 
approximation to the weak non-linear coupling scenario — 
we are able to build the exact likelihood on pixel space. 
Working in pixel space allows one to include easily non-ideal 
observational conditions, like incomplete sky coverage and 
anisotropic noise. The extension of the present work is to 
account for higher order moments in the expansion, in par- 
ticular, we are able to directly obtain constraints on gNL, 



that is the coupling parameter governing the cubic term of 
the weak non-linear expansion. 

As far as we are aware, direct constraints on gNL have 



been ma de available only very recently (Desjacques & Sel- 
jak|2009 ) and are obtained from LSS analyses: —3.5 x lO'' < 
gNL < 8.2 X 10^ at 95% CL. This constraint was obtained 
for the specific case in which the coupling parameter gov- 
erning the quadratic term (fNL) is negligible (i.e. fNL = 0, 
which is the situation required for some curvaton inflation- 
ary models, e.g. [Sasaki et al.]|2006| [Enqvist fc Takahashi| 
2008 Huang|[2009 1 . We present in this work the first direct 
measurement of gNL obtained from CMB data. In addition 
to study the particular case of ft-jL = 0, we also consider a 
more general case in which a joint estimation of fwL and gNL 
is performed. Finally, and justified by recent findings (e.g., 
Hoftuft et al.|[2009 l, we compute the N-pdf in two different 
hemispheres, and derive from it constraints on fNL and gNL 
for this hemispherical division of the celestial sphere. 

The paper is organized as follows. In Section [2] we de- 
scribe the physical model based on the local expansion of 
the CMB fluctuations and derive the full posterior probabil- 
ity, recalling how it defaults to the case already addressed 
in |Vielva fc Sanz| ( |2009[ ). In Section |3] we check the method- 
ology against WMAP-like simulations. Results on WMAP 
5-year data are presented in Section|4] Conclusions are given 
in Section [5] Finally, in Appendix |A] we provide a detailed 
computation of the full N-pdf. 



2 THE NON-GAUSSIAN MODEL 

Although current CMB measurements are well described by 
random Gaussian anisotropics (as it is predicted by the stan- 
dard inflationary scenario), observations also allows for small 
departures from Gaussianity, that could indicate the pres- 
ence of an underlying physical process generated in non- 
standard models. 



As we did in Vielva & Sanz ( 2009 1 we adopt a para- 



metric non-Gaussian model for the CMB anisotropics, that 
accounts for a small and local (i.e. point-to-point) pertur- 
bation of the CMB temperature fluctuations, around its in- 
trinsic Gaussian distribution: 



AT, = (ATOg + a [(AT,)g - ((AT, 
b[AT,)l + 0{{AT,)%). 



+ 



(1) 



The linear term ((ATi)^) is given by a Gaussian N-point 
probability density function (N-pdf) that is easily described 
in terms of the standard inflationary model. The second and 
third terms on the right-hand side are the quadratic and the 
cubic perturbation terms, respectively, and they are gov- 
erned by the a and b parameters. 

The sub-index i refers to a given direction in the sky 
that, in practice, is described in terms of a certain pixeliza- 
tion of the sphere. The operator (•) indicates averaging over 
all the pixels deflning the sky coverage. We have not con- 
sidered explicitly an instrumental noise-like term, since, for 
the particular case that we intend to explore (i.e., large-scale 
CMB data), its contribution to the measured signal (for ex- 
periments like WMAP or Planck) is negligible. Precisely at 
the large-scale regime, the term AT^ is mostly dominated 
by the Sachs- Wolfe contribution to the CMB fluctuations. 
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and can be related to the primordial gravitational potential 
$ (e.g. [Komatsu et al.||2001[ ) by: 



(2) 



where To is the CMB temperature. Small departures from 
Gaussianity of the $ potential are usually described via the 
weak non-linear coupling model: 

= $L,,+fNL - (<E-i;',))+gNL$L' +0 . (3) 

Taking into account equations [l] [2] and [3] and always con- 
sidering the specific case for scales larger than the horizon 
scale at the recombination time (i.e. above the degree scale), 
it is trivial to establish the following relations: 



fN 



To 



To 



(4) 



At this point, it is worth mentioning that the model in equa- 
tion [l] does not pretend to incorporate all the gravitational 
and non-gravitational effects, due to the evolution of the ini- 
tial quadratic potential model, but rather allows for a use- 
ful parametrization for describing a small departure from 
Gaussianity. The relationships in equation |4] have to be un- 
derstood just as an asymptotic equivalence for large scales. 

Let us simplify the notation by transforming the Gaus- 
sian term {ATi)^ into a zero mean and unit variance ran- 
dom variable (f)i . It is easy to show that equation [l] trans- 
forms into: 



+ aa I 



where: 



-AT 



1) + ba^fht + O I 



(AT), 



(5) 



(6) 



and = ((ATi)^) is the standard deviation of the CMB 
fluctuations. Trivially, the normalized Gaussian variable cp 
satisfies: 



/2n + l 



> = 





(2m - 



1)!! 



(7) 



where n ^ and m > are integer numbers, and ^ij rep- 
resents the normalized correlation between pixels i and j. 
Obviously, the N-pdf of the </> — {<j>i, (j)2, (pn} random 
field (where A'^ refers to the number of pixels on the sphere 
that are actually observed) is given by a multivariate Gaus- 
sian distribution: 



(27r)'v/2(det^)i/2 



(8) 



where ^ denotes the correlation matrix and operator de- 
notes standard matrix/vector transpose. 



As it was the case in Vielva & Sanz (2009 1, the objective 



is to obtain the N-pdf associated to the non-Gaussian x — 
{xi, X2, a;]v} field, as a function of the non-linear coupling 
parameters (or, equivalently, the a and b coefficients): 



p{x\a,b) =p{4> = 4>{x))Z. 



(9) 



In this expression, Z is the determinant of the Jacobian for 
the (p — ^ ^ transformation. Because the proposed model 



is local (i.e. point-to-point), the Jacobian matrix is diagonal 
and, therefore, Z is given by: 

d(j)i 



det 



n 



dxi 



(10) 



where: 

m = 

= 

A, = 



Both, equations [9] and [To] require the inversion of equa- 
tion [5] i.e., to express (pi as a function of Xi. After some 
algebra, it can be proved that: 

(pi ^ Xi+ riia + Vio'^ + AiCT^ + /iiO"'' + O (cr^) , (11) 

-a [xi - 1) 

2a^ — 6) — 2a^Xi 

5ab ~ 5a'^) a;* + (6a^ — Safe) x^ — a'^ 

- 2lab + 36^) x^ + (-20a* + 20a^fe) xl 

+ {&a^ -■iab)xi. (12) 

Instead of dealing with p{x\a,b), it is equivalent, but more 
convenient, to work with the log-likelihood C{x\a,b): 

A detailed computation of £(a;|a,6) is given in Appendix 
A. Let us just recall here its final expression: 

^C{x\a,b) = Fo- + (2a^ - 36 + G) a V Hcr^ 

12a''-36a'6+y6'+/^a^ (14) 

where A'' is the number of data points, and F, G, H and / are 
functions of a and b (see [AT3I ). The desired N-pdf, p{x\a,b) 
is obtained by the inversion of equation [13] and taking into 
account that p(a;|0) = p{(p = x), i.e., the known Gaussian 
N-pdf in equation [8] 



3 APPLICATION TO WMAP SIMULATIONS 

In this Section we aim to investigate the performance of the 
parameters estimation from the N-pdf derived in the previ- 
ous Section. We explore different non-Gaussian scenarios; in 
particular, we study three particular cases of special interest: 

• Case i) a 0, b — 0. This scenario would correspond, 
for example, to the case for the slow-roll standard inflation 

• Case li) a = 0, b ^ 0. This scenario would correspond 
to the particular situation for some curvaton models. 

• Case in) a 7^ 0, 6 ^ 0. It is a generic case, not repre- 
senting any specific infiationary model, but rather a general 



In particular, we will study how the determination of 
the parameters governing the non-Gaussian terms is per- 
formed, and what is the impact when one is exploring dif- 
ferent configurations. In the next subsections, we will focus, 
first, in the case in which a slow-roll standard like scenario is 
assumed (i.e., we only try to adjust for the quadratic term, 
assuming the cubic one is negligible), whatever the data is 
actually a pure quadratic or cubic model, or a general non- 
Gaussian scenario. Second, we will follow a similar analy- 
sis, but assuming the estimation of a pure cubic term. Fi- 
nally, we will address the case for a joint estimation of both 
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Figure 1. Mask at NSIDE=32 HEALPix resolution used in this 
work. It corresponds to the WMAP KQ75 mask, although the 
point source masking has not been considered, since the point 
like-emission due to extragalactic sources is negligible at the con- 
sidered resolution. At this pixel resolution, the mask keeps around 
69% of the sky. 



[quadratic and cubic) terms. In the following, we will refer 
all our results in terms of the non-linear coupling parameters 
(fNL and gNL), rather than to the a and b coefficients. 

In order to carry out this analysis, we have gener- 
ated Gaussian CMB simulations coherent with the model 
induced from t he WMAP 5-year data at NSIDE=32 
HEALPix (Gorski et al.|2005[ ) resolution (-2°). 

The procedure to generate a CMB Gaussian simula- 
tion — (Ar)g, in equation [T] — is as follows. First, we simu- 
late WMAP observations for the Ql, Q2, VI, V2, Wl, W2, 
W3, W4 difference assemblies at NSIDE=512 HEALPix res- 
olution. The Ci obtained with the cosmological parameters 
provided by the best-fit to WMAP data alone (Table 6 in 
Hinshaw et al.|2009| |, are assumed. 

Second, a single co-added CMB map is computed after- 
wards through a noise-weighted linear combination of the 
eight maps (from Ql to W4). The weights used in this lin- 
ear combination are proportional to the inverse mean noise 
variance provided by the WMAP team. They are indepen- 
dent on the position (i.e., they are uniform across the sky 
for a given difference assembly) and they are forced to be 
normalized to unity. Notice that we have not added Gaus- 
sian white noise to the different difference assembly maps, 
since we have already checked that instrumental noise plays 
a negligible role at the angular resolution in which we are 
interested (« 2°, see Vielva fc Sanz|[2009 for details). 

512 is degraded 
and a mask rep- 



Third, the co-added map at NSIDE= 
down to the final resolution of NSIDE=32, 
resenting a sky coverage like the one allowed by the WMAP 
KQ75 mask ( [Gold et al.|2009[ ) is adopted. At NSIDE=32 the 
mask keeps around 69% of the sky. The mask is given in fig- 
ure [T] Let us remark that observational constraints from an 
incomplete sky coverage can be easily taken into account by 
the local non-Gaussian model proposed in this work, since 
it is naturally defined in pixel space. This is not the case 
for other common estimators like the bispectrum, where the 
presence of an incomplete sky coverage is usually translated 
into a loss of efficiency. 

We have generated 500000 simulations of (Ar)^,, com- 
puted as described above, to estimate the correlation matrix 
$, accounting for the Gaussian CMB cross-correlations. We 
have checked that this large number of simulations is enough 



to obtain an accurate description of the CMB Gaussian fluc- 
tuations. 

In addition, we have generated another set of 1000 sim- 
ulations. These are required to carry out the statistical anal- 
ysis to check the performance of the parameter estimation. 
Each one of these 1000 (AT%j simulations are transformed 
into X (following equations [l] and [6| to study the response 
of the statistical tools as a function of the non-linear pa- 
rameters defining the local non-Gaussian model proposed in 
equation [5] 

Finally, let us remark that hereinafter the likelihood 
maximization is simply performed by exploring a grid of 
values in the parameter space of the non-linear coupling pa- 
rameters. The step used in the grid is small enough to guar- 
antee a good estimation both of the likelihood peak and 
tails. 

3.1 The recovery of Tnl in the presence of a cubic 
term 

The results obtained from the 1000 simulations are given in 
figure[2] We have explored 16 different non-Gaussian models, 
accounting for all the possible combinations obtained with 
simulated gNL values of 0, 3 x 10^, 5 x 10^ and 10®, and Inl 
values of 0, 200, 400 and 600. For each panel, we present 
the histogram of the maximum-likelihood estimation of the 
non-linear coupling quadratic parameter: fwL. Notice that 
we refer to a simulated value of a given non-linear coupling 
parameter (xnl), as xnl, whereas that its estimation via the 
maximum- likelihood is denoted as xnl . Vertical dashed lines 
in each panel, indicate the value of the maximum-likelihood 
estimation for the parameter. 

As it can be noticed from the figure, when the simulated 
data satisfies the condition of the particular explored model 
(i.e., gNL = 0, first column), the fNL is accurately and effi- 
ciently estimated, at least for values of Jnl < 600. Actually, 
this is a result that we already obtained in |Vielva fc Sanz| 
(20091, which indicates that for Jnl > 600, the perturbative 



model stars to be not valid any longer. 

However, when the simulated non-Gaussian maps also 
contain a significant contribution from a cubic term, the 
bias in the determination of the fNL parameter stars to be 
evident already for lower values of the simulated Inl. It is 
interesting to notice that, even if the simulated gNL is large 
(for instance gNL = 10®), we do not see any significant bias 
in Inl, for simulated Inl values lower than 200. 

Summarizing, we can infer that for non-Gaussian sce- 
narios with |fNL| 400 and |gNL| J; 5 x 10'', no significant 
bias on the estimation of a pure quadratic term is found. It 
is worth mentioning that these range of values are in agree- 
ment with predictions from most of the physically motivated 
non-Gaussian infiationary models. Notice that, in general, 
even for the cases in which a bias is observed, the efficiency 
in the determination of fNL (somehow related to the width 
of the histograms) is almost unaltered. 

3.2 The recovery of gNL in the presence of a 
quadratic term 

As for the previous case, a graphical representation of the 
results obtained from the 1000 simulations is given in fig- 
ure [3] We have explored the same 16 different non-Gaussian 
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Figure 2. These panels represent the accuracy and efficiency on the estimation of the f^L parameter. From left to right, the columns 
correspond to simulated g^L values of: 0, 3 X 10^, 5 X 10^ and 10®. Similarly, from top to bottom, rows correspond to simulated f^jL 
values of: 0, 200, 400 and 600. The histograms show the distribution of the obtained values of ff,jL for each case. The vertical dashed 
lines indicate the simulated fjsjL value, and help to identify the presence of a possible bias. 
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Figure 3. These panels represent the accuracy and efficiency on the estimation of the gjML parameter. From left to right, the columns 
correspond to simulated g^L values of: 0, 3 X 10^, 5 X 10^ and 10®. Similarly, from top to bottom, rows correspond to simulated f^jL 
values of: 0, 200, 400 and 600. The histograms show the distribution of the obtained values of gujL for each case. The vertical dashed 
lines indicate the simulated g^L value, and help to identify the presence of a possible bias. 
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models already described above. As it can be noticed from 
the figme, when the simulated data corresponds to the ex- 
plored model (i.e., Tnl = 0, first row), the gNL parameter is 
reasonably estimated, at least for simulated gNL < 10®. 

However, when the simulated non-Gaussian maps also 
contain a significant contribution from the quadratic term, 
a bias in the determination of the gNL parameter stars to be 
notorious for lower values of the simulated gNL coefficient. 
In particular, the plots of the first column (i.e., gNL = 0) 
show a clear bias on gNL. This indicates that, when the an- 
alyzed case corresponds to a pure quadratic scenario, and a 
pure cubic model is assumed, the gNL estimator is sensitive 
to the quadratic non-Gaussianity and, somehow, it absorbs 
the non-Gaussianity in the form of a fake cubic term. In 
particular, an input value of Tnl = 400 is determined as a 
pure gNL ~ 2.5 x 10^. Notice that this was not the situ- 
ation for the previous case, where the fNL estimation was 
not sensitive to the presence of a pure cubic model (at least 
for reasonable values of gNL). This is an expected results, 
since, any skewned distribution would imply the presence 
of a certain degree of kurtosis, whereas the opposite is not 
necessary true. 



to the inverse average noise variance across the sky, and 
are normalized to unity. This linear combination is made 
at NSIDE=512 HEALPix resolution, being degraded after- 
wards down to NSIDE=32. 

Under these circumstances, we are in the same condition 
as for the analysis performed on the simulations described 
in Section |3] Therefore, the theoretical multinormal covari- 
ance of the CMB temperature fluctuations (^) is the one 
already computed with the 500000 simulations (see previ- 
ous Section). 

Two different analysis were performed. The first ac- 
counts for an all-sky study (except for the sky regions cov- 
ered by the Galactic mask described in the previous Sec- 
tion) , where constraints on the non-linear coupling parame- 
ters from different scenarios are presented. We will present as 
well results derived from a model selection approach, where 
we investigate which are the models that are more favoured 
by the data. The second analysis explores the WMAP data 
compatibility with the local non-Gaussian model in two dif- 
ferent hemispheres. In particular, we have studied indepen- 
dently the two hemispheres related to the dipolar pattern 



described in Hoftuft et al. ( 2009 1 



3.3 The general case: the joint recovery of fNL and 

gNL 

Finally, we have also explored the case of a joint estima- 
tion of the quadratic and cubic terms. The results obtained 
from the 1000 simulations are given in figure |4] As for the 
previous cases, we have explored the same 16 different non- 
Gaussian models already described above. The plots repre- 
sent the contours of the 2-D histograms obtained for the pair 

^fNLjgNL^ of the maximum-likelihood estimation. Vertical 

and horizontal dashed lines indicate the simulated Jnl and 
gNL values, respectively. 

As it can be noticed from the figure, only for the regime 
|fNL| 400 and |gNL| 5 x lO'', we obtain an accurate 
and efficient estimation of the non-linear coupling param- 
eters. As it was reported above, this regime correspond to 
the boundaries obtained from the pure fNL case. 

It is interesting to notice the presence of very large bi- 
ases for cases outside of the previous range. In particular, 
estimations tend to move towards a region of the parameter 
space of larger values of both, fNL and gNL . Only a secondary 
peak in the 2-D histogram corresponds to the simulated pair 
of values. 

This result, combined with the previous ones, indicates 
that the non-Gaussian model proposed in equation [l] is only 
valid up to values of the quadratic and cubic terms of around 
1% and 0.05%, respectively. 



4 APPLICATION TO WMAP 5- YEAR DATA 

We have studied the compatibility of the WMAP 5-year data 
with a non-Gaussian model as the one described in equa- 
tion [T] In particular, we have analyzed the co-added CMB 
map generated from the global noise-weighted linear com- 
bination of the reduced foreground maps for the Ql, Q2, 
VI, V2, Wl, W 2, W3 and W4 diflFerence assemblies (see 
Gold et al.||2009[ for details). The weights are proportional 



4.1 All-sky analysis 

We have computed the full N-pdf in equation[9] for three dif- 
ferent scenarios: a non-Gaussian model with a pure quadratic 
term (i.e., gNL = 0), another case with a pure cubic term 
(i.e., fNL = 0), and a general non-Gaussian model (i.e., 
fNL / and gNL / 0), 

Results are given in figure [5] Left panel shows the like- 
lihood obtained for the first case: p(a;jfNL). Actually, this 
result is the one that we already obtained in our previous 
work ( jVielva fc Sanz||2009| |. The maximum-likelihood esti- 
mation for the quadratic factor is Inl = The con- 
straint on the non-Gaussian parameter is: —154 < Inl < 94 
at 95%. 

The middle panel in figure [5] presents the likelihood 
obtained from a model with a pure cubic term: p(a;|gNL). 
The maximum-likelihood estimation for the quadratic fac- 
tor is gNL ~ 42785. The constraint on the parameter is: 
-5.6 X 10^ < gNL < 6.4 X 10* at 95%. This result is com- 
patible with a previous finding obtained from the analysis of 



LSS data (Desjacques & Seljak 20091. The result reported 



in this work is, as far as we know, the first direct constraint 
of gNL from CMB data alone. 

The right panel in figure [5] shows the contour levels 
at the 68%, 95% and 99% CL, for the likelihood obtained 
from an analysis of a general quadratic and cubic model: 
p (a; I fNL, gNL). Notice that the maximum likelihood estima- 
tion for the fNL and gNL parameters are similar to those 
obtained from the previous cases (where the pure models 
were investigated). Even more, the marginalized distribu- 
tions for the two parameters are extremely similar to the 
likelihood distributions discussed previously, and, therefore. 



Notice that in [Vielva fc Sanz] | |2009| we used a different defi- 
nition between the primordial gravitational potential $ and the 
CMB temperature. This difference implies that, in our previous 
work, the Tnl parameter had an opposite sign with respect to the 
definition used in this paper. 
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Figure 4. These panels represent the accuracy and efficiency on the joint estimation of the ff,jL and g^L parameters. From left to right, 
the columns correspond to simulated gNL values of: 0, 3 X 10^, 5 X 10^ and 10®. Similarly, from top to bottom, rows correspond to 
simulated Inl values of: 0, 200, 400 and 600. The contours show the distribution of the obtained values of the pair ^fNLigNL^ for each 

case. The vertical and horizontal dashed lines indicate the simulated Inl a-nd g^L values, respectively, and help to identify the presence 
of a possible bias. 
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Figure 5. Likelihood distribution of tlie non-linear parameters obtained by analyzing 1000 simulations, according to the local non- 
Gaussian model given in equation[l] Left panel correspond to a pure fjML analysis: p (£c|fNL)- Middle plot shows the result for a pure gNL 
model: p(a'|gNL)- Right panel provides the 68%, 95% and 99% contour levels of the likelihood obtained from a joint fjMLi gNL analysis 
of the WMAP 5-year data: p (a;|fNL, Snl)- 



the constraints on the non-linear coupling parameters are 
virtually the same. 

Finally, we want to comment a few words about two 
issues: the incorporation of possible a prion information re- 
lated to the parameters defining the non-Gaussian model, 
and the application of model selection criteria (or hypoth- 
esis tests) to discriminate among the Gaussian model and 
different non-Gaussian models. 

As we largely discussed in our previous work (Vielva fcj 
Sanz||2009[ ), one of the major advantages of computing the 
full N-pdf on the non-Gaussian model is that, in addition 
to provide a maximum-likelihood estimation for the non- 
linear coupling parameters, we have a full description of the 
statistical properties of the problem. More in particular, if 
we could have any physical (or empirical) motivated prior 
for the fNL and gNL parameters, it could be used together 
with the likelihood function to perform a full Bayesian pa- 
rameter estimation. This aspect has not been considered in 
this work, precisely because such a well motivated prior is 
lacking. Actually, a possible and trivial a priori information 
that could be used in this specific case, would be to limit 
the range of values that can be taken by fNL and gNL, such 
as the non-Gaussian model is, indeed, a local perturbation 
of a Gaussian field (i.e., the typical values that we discussed 
in Section[3|. However, these priors do not seem to be quite 
useful since, first, we do not have any evidence to chose any 
different form for the prior that an uniform value over the 
parameters range; and, second, the limits of these ranges 
are somehow arbitrary. These kind of priors do not provide 
any further knowledge on the Bayesian parameter determi- 
nation: as it is well known, such estimation would be totally 
driven by the likelihood itself, since it is fully defined within 
any reasonable a priori ranges. 

The possibility of performing a model selection ap- 
proach is an extra advantage of dealing with the full N-pdf. 
Of course, under the presence of an hypotetical well moti- 
vated prior on the non-linear coupling parameters, model 
selection could be done in terms of the Bayesian evidence 



or the ratio of posterior probabilities (see Vielva & Sanz 



2009 for a specific discussion on this application) . However, 



the lack of such a prior (as we discussed above), makes the 
application of a full Bayesian model selection framework sig- 
nificantly less powerful than in other situations: as it is very 
well known, the use of uniform priors for all the parameters 
would provide very little information, since the results would 



be very much dependent on the size of the parameters range. 
Despite this, we can still make a worthy use of the likelihood 
to perform model selection. In particular, some asymptotic 
model selection criteria, like the Akaike Information Crite- 
ria (AIC, [Akaike||1973[) an d the Bayesian Information Cri- 
teria (BIG, Schwarz 19781, can be applied. Both methods 
provide a ranging index for competitive hypotheses, where 
the most likely one is indicated by the lowest value of the 
index. The AIC and BIG indices depend on the maximum 
value of the log-likelihood (max [C (a;|0)] = £■): 



AIC (HO = 2lp-C 



BIG (HO = 2(^1 log AT 



where p is the number of parameters that determine the hy- 
pothesis or model Hi. We have applied these two asymptotic 
model selection criteria to the WMAP 5-year data. Defin- 
ing the Gaussian model as Ho, the pure quadratic model 
as Hi, the pure cubic model as H2, and the general non- 
Gaussian model as H3 , and considering the maximum value 
for the log-likelihoods obtained for all these cases, we ob- 
tain: AIG(//o) < AIG(H"i) < AIC (Hi) < AlCiHi), and 
BIG (Ho) < BIG (Hi) < BIG {H2) < BIG {H3). That is, the 
most likely model is the Gaussian one (what is in agreement 
with the results obtained from the parameter determination, 
since fNL = and gNL = can not be rejected at any mean- 
ingful confidence level). Among the non-Gaussian models, a 
pure fNL model is the most likely scenario, being a joint fNL, 
gNL model the most disfavoured by the WMAP 5-year data. 

4.2 Hemispherical analysis 

Among the large number of the WMAP anomalies that have 
been reported in the literature, an anisotropy manifested in 
the form of a hemispherical asymmetry is one of the topics 



that has been more largely studied (e.g., Eriksen et al. 2004a 



Hansen et al. 2004b I. Most of the works related to this is- 



sue, have reported that such asymmetry is more marked for 
a north-south hemispherical division relatively close to the 
Northern and Southern Ecliptic hemispheres. 

In a recent work, Hoft uft et al.| ( |2009[ ) reported that 
large scale WMAP data was compatible with such kind of 
anisotropy, in the form of a dipolar modulation defined by a 
preferred direction pointing toward the Galactic coordinates 
{l,b) = (224°, -22°). 
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Figure 6. These plots show the areas of the sky that are independently analyzed. The panel on the left accounts for the sky that, 
being allowed by the Galactic mask (see Section [sj, corresponds to the Northern hemisphere of the sky division considered in this work. 
Equivalently, the right panel presents the region of the sky that is analyzed when the Southern hemisphere is addressed. 



Motivated by these results, we have repeated the anal- 
ysis described in the previous subsection, but addressing in- 
dependently the two hemispheres associated to the dipolar 



pattern found by Hoftuft et al. (20091. Hereinafter, we will 



refer to the Northern hemisphere of this dipolar pattern, as 
the half the of celestial sphere whose pole is closer to the 
Northern Ecliptic Pole, and, equivalently, we will indicate 
as the Southern hemisphere the complementary half of the 
sky. The corresponding areas of the sky that are analyzed are 
shown in figure[6] The left and right panels show the allowed 
sky regions, when the Northern and Southern hemispheres 
of the dipolar pattern are independently addressed, respec- 
tively. Notice that the regions not allowed by the Galactic 
mask are also excluded from the analysis. The portion of 
the sky that is analyzed is around 34% for the Northern 
hemisphere, and around 35% for the Southern half. 

As we discussed in the previous subsection, the con- 
straints of the fNL and gNL parameters obtained from the 
analysis of pure quadratic and cubic non-Gaussian models, 
do not differ significantly from those obtained from a general 
analysis of a joint scenario. This is expected for a regime of 
relatively low values of the non-linear coupling parameters. 
For that reason, in the present study we will only consider 
the following two cases: a pure quadratic (i.e. gNL = 0) and 
a pure cubic (i.e. fwL = 0) models. Results are given in fig- 
ure [7] We present the likelihood probabilities for the first 
case — p (a; I fNL) — in the left plot, and the one correspond- 
ing to the second case — p (a;|gNL) — in the right panel. Each 
plot shows the results for the Northern (dashed lines), and 
the Southern (dot-dashed lines) hemispheres. The maximum 
likelihood estimation for the non-linear coupling parameters 
are given as vertical lines. 

The right panel shows that both hemisphere have a sim- 
ilar likelihood (p (a;|gNL)) for the case of a pure cubic model. 
However, interestingly, it is not the case when addressing a 
gNL = scenario. In this case, we notice two important re- 
sults. First, whereas the fNL estimation from the analysis of 
the Northern hemisphere provides a constraint compatible 
with the Gaussian scenario, it is not the case for the South- 
ern hemisphere. In fact, we find that fNL < at 96% CL. In 
particular we find: fNL ~ —164 ± 62. Second, the distance 
between both distributions is too large. Let us make use of 
the KuUback-Leibler divergence (KLD, [Kullback fc Leibler| 



Dn,s = / dfNLPn (a;jfNL)log 



(15) 



likelihoods p„ (ccjfNL) and ps (a;jfNL): 

Pn (a; I fNL) 

Ps (a;|fNL) ' 

where p„ (a;|fNL) and ps (a;|fNL) are the likelihoods for the 
Northern and the Southern hemispheres, respectively. Actu- 
ally, we use the symmetrized statistitic D, defined as: 



(16) 



We have found that the distance D for the likelihood dis- 
tributions of the fNL parameter estimated in the Northern 
and the Southern hemispheres defined by the dipolar pat- 



tern described by Hoftuft et al. ( 2009 1 is much larger than 



it would be expected from Gaussian and isotropic random 
CMB simulations. In particular, such a distance has a p- 
value « 0.04. This result is a further evidence on the largely 
discussed WMAP North-South asymmetry, and it is as well 
an indication that such asymmetry is manifested in terms of 
the non-Gaussianity of the CMB temperature fiuctuations, 



in agreement with previous results (e.g.. Park 
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1951 1 as a measurement of the distance between the two 



At this point, it is worth recalling that the analysis de- 
scribed above has been performed assuming isotropy, i.e., 
we have used the same type of correlations to described the 
second-order statistics in both the Northen and in the South- 
ern hemispheres. However, the result obtained by |Hoftuft et| 
al. (2009j) indicates that these two hemisphere might be de- 
scribed by two different correlations (i.e., the sky would not 
be isotropic any longer). The dipolar modulation proponed 
by Hoftuft et al. ( 2009 1 was small (its amplitude was lower 
than 0.7%), but significant (a 3.3(t detection was claimed). 
Assuming this point, we have repeated our previous analy- 
sis, but using different statistical properties for the correla- 
tion matrices in the two halves. The way we have estimated 
these new correlation matrices is as follows: we have gen- 
erated 500,000 simulations (in the same way it has been 
already described at the beginning of Section [4] and, once 
the co-added maps are degraded to NSIDE=32, each one of 
the simulations have been modified by applying the dipo- 
lar modulation estimated by [Hoftuft et al.| (2009) from the 
WMAP data. It is from these modulated simulations that 
we have estimated the new correlation matrices needed to 
estimate the likelihood probabilities. The result of this test 
is presented in figure [S] The conclusions related to the fwL 
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Figure 7. The left panel present the likelihood on fj^L obtained from a pure quadratic analysis — pi^l^Nh) — i whereas the right plot 
provides the likelihood on g^L from a pure cubic study — ^p(3;|gnjL)- Dashed lines correspond to the Northern hemisphere, whereas dot- 
dashed lines are for the Southern half. Vertical lines indicate the maximum likelihood estimation of the non-linear coupling parameters: 
fNL and gNL- 



estimation are essentially the same: on the one hand the 
analysis of the Northern hemisphere provides a constraint 
compatible with the Gaussian scenario, whereas the South- 
em hemisphere is not; on the other hand, the distance be- 
tween both distributions is again too large (it also has a p- 
value -as compared to, in this case, anisotropic simulations- 
increases up to « 0.09 (estimated in terms of the KLD). 
However, despite this slight change in the fNL hemispheri- 
cal estimation, dramatic differences can be observed for the 
pure cubic scenario. Interestingly, accounting for the dipo- 
lar modulation correction reveals an extra departure form 
anisotropy related to the gNL constraints. The dipolar mod- 
ulation makes the maximum likelihood estimation of gNL 
highly incompatible between both hemispheres. In partic- 
ular, the distance between both distributions is extremely 
rare as compared with the expected behaviour from Gaus- 
sian and anisotropic CMB simulations (generated, as ex- 
plained above, by applying the dipolar modulation reported 
by |Hoftuft et al.| ( |2009[ )): it has a p-value of « 0.002, in the 
sense of the KLD. 

We have also studied whether a WMAP data corrected 



by the dipolar modulation found by Hoftuft et al. ( 2009 1 



could present a behaviour compatible with the Gaussian and 
isotropic hypotheses. Results for the corrected WMAP data 
are given in figure |9] Notice that we do not see any signif- 
icant differences from the previous situation (i.e., the case 
in which uncorrected WMAP data was analyzed assuming 
anisotropy): the dipolar modulation does not affect to the 
fNL and gNL constraints. 

Summarizing, the results obtained in this subsection 
seem to confirm that there is some kind of anomaly related to 
an hemispherical asymmetry as the one defined by the dipo- 



lar pattern reported by Hoftuft et al. (20091, in the sense of 



the fNL parameter. Even more, when WMAP data is ana- 
lyzed using correlations compatible with the dipolar modula- 



tion suggested by Hoftuft et al. ( 2009 1, not only asymmetries 
related to the fNL parameter are clear, but also associated 
to the cubic term (i.e., the gNL parameter). Intriguingly, the 
correction of the WMAP data in terms of this dipolar mod- 
ulation is not enough to obtain a CMB signal compatible 
with a Gaussian and isotropic random field. 



At this point, it is worth mentioning that the dipolar 
modulation of Hoftuft et al.| (12009 ) was obtained by consid- 
ering second-order moments of the CMB data and, there- 
fore, this correction only addresses the problems related to 
an asymmetry in terms of this order. Hence, it is not to- 
tally surprising that this dipolar modulation correction is 
not sufficiently satisfactory to solve the anomaly reported 
in this work, since such anomaly is related to higher-order 
moments. It is also need to point out that we have tested 
that the dipolar modulation correction of the WMAP data 
does not affect the results obtained from an all-sky analysis 
of the CMB data. 

Finally, let us recall that, in a recent work, [Rudjord et] 
|al.| ( |2009[ ) searched for specific asymmetries related to the lo- 
cal estimation of the fNL parameter, by using needlets. Con- 
trarily to our findings, in this work no significant asymmetry 
was found when analyzing WMAP data. There are some dif- 
ferences between the analyses that could explain the discrep- 
ancy, although they have to be taken as mere suggestions. 
First, the kind of non-Gaussianity that is probed by each 
work is different: whereas the Rudjord et al. (2009) paper ex- 
plore a fNL model that is local in the gravitational potential 
(from which the non-Gaussian temperature fluctuations are 
obtained taking into account all the gravitational effects), 
here we adopt a local model in the Sachs- Wolfe regime. Sec- 
ond, they work at the best WMAP resolution (around 10-20 
arcmins), whereas we focus on scales of around 2°. Third, 
we explore an specific division of the sky (the one reported 
by [Hoftuft et al.] ( |200 9)), whereas they consider several di- 
visions that, not necessarily, have to match the one used by 
us (they explore hemispherical divisions within an interval 
of around 30°). 



5 CONCLUSIONS 



We have presented an extension of our previous work ( |Vielva| 
&: Sanz|2009 1, by defining a parametric non-Gaussian model 
for the CMB temperature fluctuations. The non-Gaussian 
model is a local perturbation of the standard CMB Gaus- 
sian field, which (under certain circumstances) is related to 
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Figure 8. As in figure [7] but for the case in which the WMAP data has been analyzed using correlation matrices that account for the 
dipolar modulation. 
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Figure 9. As in figure [7] but for the case in which the WMAP data has been corrected by the dipolar modulation. 



an approximative form of the weak non-linear coupling in- 
flationary model at scales larger than the horizon scale at 
the recombination time (i.e., above the degree scale, see for 
instance Komatsu et al.|f2001 Liguori et al.|[2003 |. For this 
model, we are able to build the posterior probability of the 
data given the non-linear parameters fNL and gNL- From 
these pdfs, optimal maximum likelihood estimators of these 
parameters can be obtained. 

We have verified with WMAP-like simulations that the 
maximum likelihood estimation of the quadratic non-linear 
coupling parameters (fNL) is unbiased, at least for a reason- 
able range of values, even when non-Gaussian simulations 
also account for a cubic term. In particular, we found that 
for simulated non-Gaussian coefficients such as I^nlI J; 400 
and I gNL I ^ 5 X 10^, the estimation of fNL is accurate and 
efficient. However, when trying to study the case in which 
only a pure cubic model is addressed, the situation is dif- 
ferent. In particular, the simulated quadratic term has an 
important impact on the estimation of gNL- For instance, 
if a pure quadratic non-Gaussian model is simulated with a 
value of Jnl = 400, a value of gNL ~ 2.5 x 10^ is wrongly es- 
timated. This results indicates, obviously, that the quadratic 
term is more important than the cubic one in the expansion 
of the local non-Gaussian model, and, therefore, that not 
accounting properly for the former might have a dramatic 



impact on the latter. Contrarily, the opposite situation is 
much more unlikely. Finally, we have investigated the joint 
estimation of the fNL and gNL parameters. In this case we 
find that for a similar regime as the one mentioned above 
(i.e., |fNL| ^ 400 and |gNL| 5 x 10^), an accurate and ef- 
ficient estimation of the non-linear coupling parameters is 
obtained. However, for larger values of these coefficients, we 
find that the parameter estimation is highly biased, favour- 
ing a region of the parameter space of larger values for both 
coefficients. 

We have addressed, afterwards, the analysis of the 
WMAP 5-year data. We have consider two different anal- 
yses. First, we have investigated the case of an all-sky anal- 
ysis (except for the Galactic area not allowed by the WMAP 
KQ75 mask). Second, and motivated by previous findings, 
we have performed a separated analysis of two hemispheres. 
In particular, the hemispherical division associated to the 



dipolar pattern found by Hoftuft et al. ( 2009 1 was consid- 
ered. 

Regarding the all-sky analysis, we find, for the case in 
which a pure quadratic model is investigated, the same result 



that we already found in our previous work (Vielva & Sanz 



20091. In particular, we determine that —154 < Inl < 94 at 
95%. Equivalently, for the case of a pure cubic non-Gaussian 
model we estabUsh -5.6x lO'^ < gNL < 6.4x 10^ at 95%. This 
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is in agreement with a recent work by |Desjacques fc Seljak| 
(20091, where an analysis of LSS data was performed. The 



result that we provide in this paper is, as far as we know, the 
first direct constraint on gNL from CMB data alone. Finally, 
we have also investigated the case of a joint estimation of the 
quadratic and cubic non-linear coupling parameters. In this 
case, the constraints obtained on Inl and gNL are virtually 
the same as the ones already reported for the independent 
analyses. 

We have performed a model selection to evaluate which 
of the four hypotheses (i.e. Gaussianity, a pure quadratic 
model, a pure cubic model and a general non-Gaussian sce- 
nario) is more likely. Since a well motivated prior for the non- 
linear coupling parameters is lacking, we have used asymp- 
totic model selection criteria (like AIC and BIG), instead of 
more powerful Bayesian approaches, like the Bayesian Evi- 
dence or the posterior ratio test. Both methodologies (AIG 
and BIG) indicates that the Gaussian hypothesis is more 
likely than any of the non-Gaussian models. We also found 
that, among the non-Gaussian scenarios, the one with a pure 
quadratic model is the most favoured one, whereas the gen- 
eral one (i.e. fNL 7^ and gNL 7^ 0) is the most unlikely. 

The analysis of the WMAP data in two hemispheres 
revealed that, whereas both halves of the sky present sim- 
ilar constraints on the gNL parameter (and, in both cases, 
not indicating any significant incompatibility with the zero 
value), the analysis of a pure quadratic scenario showed a 
clear asymmetry. First, the Inl value in the hemisphere 
whose pole is closer to the Southern Ecliptic Pole, is fNL = 
— 164± 62. Which implies that fNL < at 96%. Even more, 
the distance between both likelihoods (given in terms of the 
KLD) presents a p- value of « 0.04. 

We have also analyzed the WMAP data after by con- 
sidering different correlation properties in each hemisphere 



(according to the dipolar modulation described by Hoftuft et 



al. ( 2009 1). We have tested that the behaviour found for the 
fNL is practically the same as before, and that, in addition, 
an extra anomaly appears associated with the gNL param- 
eter. In particular, the distance between both likelihoods 
is anomalously large as well (it corresponds to a p-value of 
< 0.002. Further test was performed after correcting WMAP 
data from the dipolar modulation. In this case the asymme- 
tries in the maximum-likelihood estimations of both non- 
linear coupling parameters remain unaltered. Hence, these 
results indicate that, as it has been previously reported in 
other works, there are evidences of anisotropy in the WMAP 
data, reflected as an asymmetry between two opposite hemi- 
spheres. Such anomaly is related to a different distribution 
for a non-linear coupling parameter related to the quadratic 
term. However, a correction in terms of a dipolar modula- 
tion as the one proposed by Hoftuft et al. (20091, seems not 



to be sufficient to account for this anomaly related to the 
likelihood distribution of the fNL parameter. 
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APPENDIX A: N-PDF DERIVATION 

In this Appendix we provide a detailed derivation of the N- 
pdf for the non-linear local model given by equation [5] As 
it has been already discussed, the objective is to make use 
of this likelihood to constrain the parameters that define 
the perturbative model. In particular, we are interested in 
estimating the parameters governing the quadratic (a) and 
cubic (fe) terms, that, as it was explained in Section [2] can 
be related (under certain circumstances) to the non-linear 
coupling inflationary parameters fNL and gNL, respectively. 

Let us recall here the expression for the N-pdf for the 
non-Gaussian model, as a function of the multinormal N-pdf 
(p (</)), equation |9|: 



p{x\a,b) =p{ct> = 4>{x))Z, 
where Z is the determinant of the Jacobian of the 



(Al) 



transformation (equation 10 1. For practical purposes (i.e., 
in order to constrain the non-linear coupling parameters) it 
is more convenient to deal with the log-likelihood £ {x\a, b) 



(given in equation 13 1, instead of p{x\a,b) (given in equa- 
tion |9|. Obviously, the later can be easily obtained by the 
inversion of equation |13| 

Replacing equation[9]into equation |13| it is easy to show 
that the computation of C {x\a, b) implies to solve two terms: 



£.{x^\a,b) =logZ-fn, 



(A2) 



^og ^^fc^of^ ' s.'^'i' recall it, p(a;10) 



where Q 

p{4> = x)- 

Let us address the determination of these two terms 
(log Z and Q.) independently. 

Al Determination of the log-Jacobian 

According to equation [To] the log-Jacobian for the (f) — ^ ^ 
transformation is given by: 



\ogZ = Y.^ 



(A3) 



where ifi as a function of Xi is given in equation |11[ It is 
straightforward to prove that: 



dxi 



= 1+2^ 9m,iO , 

m — 1 



(A4) 



where 



~2axi 

<?2,j = 3 (2a^ — b) xj — 2a^ 

53, » = 20 {-a^ ^ ab) 4 + 6 (2a^ - ab) Xi 

g^^i = 5 (I4a^ - 21a^6 + 36^) xt + 60(-a'' + a'^b)xi 



+2, (2a'' - a^b) 



(A5) 



Since we are considering a perturbative non-Gaussian 
model, then X]m=i 9m,icr"^ <^ 1 and, therefore, the log- 
Jacobian can be easily derived from the Taylor expansion 
for the logarithm function: 



logZ ^igi,i)a 



92,. 



1 2 \ 2 



(A6) 



+ 33 



1 3 



51,^52, 



+ 34 



i^32,i + 5i,«52,i - fli,ifl3,i - ].g\,., 



Taking into account that: 
1 ^ 

— ^ = 1 -f 2 (a" + 3&) + C» (a^) 
1=1 

1 ^ 

-^x? = 6aa + C(a=') 



JV 



N ^ 



3 + Ola' 



(A7) 



(as it can be trivially showed from equations [5] and [tJi and 
replacing equation | A5| into equation |A6[ one finally gets: 



j^logZ^ (2a^ - 36) cr^ + ^12a* - 36a^6 + y ) (A8) 
up to the appropriate order in a. 
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A2 Determination of Q, 

Since we are dealing with a perturbative model (i.e., <p — 
X + e, with e <C a; and where e = rfa + v>a^ + /itj'^ + Act*, 
according to equation 111 it is interesting to notice that the 
probability function p{cf> = cf) (x)) can be expanded in terms 
of the Taylor expansion, up to the appropriate order: 

Pi4> = 4>{x)) = P + e'd,P+^e'e'd,,P+^e'e'e''d,ikP 

+ ^6'e^e^'9,,MP + 0(5), (A9) 

where, for simplicity, we define P=p{(f> = x + 0)=p {x\0). 
The function diY (x) is the derivative of Y (x) with re- 
spect to Xi, i.e.: diY (x) = Equivalently, dijY (x) = 



8: 
dx 



^d,Y{x), d,,kY{x) = ^k,Y{x), and dijkiY (x) 



^ ^^JkY (x) 



Taking into account this expansion, we have:: 



p{x\0) 2 6 

+ ^A'''''e'e'e''e\ (AlO) 

where. A' = A'' = A'^'' = and A'-'" = 

^'^p' "'^ . It is trivial to show that these quantities are related 
to the data x, its expected correlation ^ — [^ij], and the 
model (i.e., a and b). After some algebra, one trivially finds: 

A, = -^-,'x\ (All) 
Aij = + AiAj, 

Aijk = AijAk — j4j^jj. — ^iCjfe 
^ijfci = AiAjAkAi + + Cifc^Cj/ + ^jk(,ii ^ 

-AkAi^-/ - A,A^-,' - A^Ai^-,' 
—AkAi^ji — AkAj(^^i — AiAj(,^.i . 

Since we are interested in the quantity log ^'•^^'^q^-'-' , the 
logarithm function can be expanded up to the appropriate 
order. It is straightforward to obtain: 

n = Fa + Ga^ + Ha^ + la^ + O (5) , (A12) 

where: 

F = -^(5^.), (A13) 



1 /c, 1 

N y^"^ + 2' 



We have defined the operator Sa/s = a^£,r/l3^ , and tj, v, A 
and iJ, where defined in equation [12] 



Therefore, taking into account equations |A8| and |A12[ 
the final expression for the log-likelihood £ {x\a,b) is: 

j^C{x\a,b) = Fa + (2q,^ - 36 + G) + ffa^ 

+ M2a'' ^ 36a^6+ yfc^ + Jj a*. (A14) 



